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1.  Summary 

We  have  developed  a  new  numerical  solver  that  models  propagation  of  femtosecond  pulses  in  multi¬ 
photon  absorbing  media.  The  Mathematica- based  numerical  calculation  predicts  how  transmission 
through  a  medium  comprising  two-photon  absorbing  (2PA)  chromophores  at  a  given  concentration 
and  with  given  molecular  parameters  changes  the  energy,  the  temporal-  and  spectral  amplitude  and  the 
phase  of  the  pulse.  Our  formalism  makes  use  of  full  three-level  stochastic  density  matrix  (SDM) 
equations  of  motion,  which  allows  us  to  study  optical  power  limiting  medium  in  a  broad  range  of 
parameters,  including  the  coherent  regime  of  the  light-matter  interaction,  i.e.  when  the  duration  of  the 
pulse  is  comparable  or  shorter  that  the  decoherence  time  of  the  medium.  Furthermore,  because  the 
SDM  equations  allow  modeling  of  arbitrary  non-Loretzian  line  shapes,  our  calculation  takes  into 
account  the  realistic  resonance  line  shape,  including  both  one-photon  and  two-photon  spectra  of  the 
chromophore. 

We  apply  the  new  numerical  solver  to  100-fs  Gaussian  pulse  of  variable  peak  intensity  10-10  GW 
cm2.  The  pulse  is  incident  on  two-level  or  three-level  medium  where  the  molecular  parameters, 
including  the  transition  frequencies,  the  values  of  the  transition  electric  dipole  moments  and  the 
difference  between  the  permanent  electric  dipole  moments  are  modeled  after  a  well-known  2PA 
chromophore  trans-4,4’-Bis(diphenylamino)stilbene  (BDPAS).  Modeling  of  most  other  2PA 
chromophores  is  readily  feasible  by  appropriate  change  of  the  parameters.  With  the  goal  of  validating 
our  calculation  we  verify  that  at  moderate  intensity  <10  GW  cm  the  numerical  results  are 
quantitatively  in  agreement  with  the  analytical  solution.  We  calculate  the  transmission  of  100-fs  pulse 
through  BDPAS  medium  at  chromophore  concentration  0.5  10  cm'  and  show  that  the  minimum 
propagation  distance  for  the  optical  power  limiting  to  take  effect  is  about  0.4  cm  and  that  larger 
distance  >0.5  cm  is  needed  for  incident  intensity  >3  10  GW  cm'  due  to  saturation  of  the  2PA.  At  the 
same  time,  the  temporal  shape  of  the  transmitted  pulse  changes  from  flat-topping  at  moderate  intensity 
to  an  asymmetrical  shape  at  high  intensity,  whereas  the  spectral  maximum  shifts  to  higher  frequency. 
These  results  agree  with  earlier  experimental  observations  in  quantum-dot-doped  waveguides.  Finally, 
we  calculate  few-cycle  5-fs  pulse  propagating  in  dipolar  two-level  medium  and  show  that  the  shape 
changes  depending  on  the  carrier-envelope  phase  (CEP)  of  the  incident  pulse. 

We  envisage  that  the  solver  can  be  applied  for  simulation  of  a  broad  range  of  physical  problems 
including  coherent-  and  incoherent  regimes  of  optical  power  limiting,  saturation,  CEP  effects, 
electromagnetically  induced  transparency,  soliton  formation  etc. 
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2.  Introduction 

Interaction  of  high  peak  intensity  light  pulses  with  a  multi-photon  absorbing  medium  is  an  important 
issue  in  a  number  of  ultrafast  photonics  applications.  Developing  a  practical  yet  accurate  numerical 
tool  for  modeling  of  the  propagation  of  a  femtosecond  laser  pulse  in  a  two-photon  absorbing  (2PA) 
media  would  benefit  optimization  of  optical  power  limiting  devices,  all-optical  switches,  two-photon 
fluorescence  microscopes  etc. 

In  this  project  we  have  developed  a  Mathematica- based  numerical  calculation  that  simulates  how  the 
energy,  temporal-  and  spectral  amplitude  and  phase  of  a  femtosecond  pulse  develops  as  the  pulse 
propagates  in  a  medium  comprising  two-photon  absorbing  chromophores.  We  use  density  matrix 
equation  formalism,  which  allows  us  to  study  the  transmission  through  an  optical  power  limiting 
medium  in  coherent  regime,  where  the  duration  of  the  pulse  is  comparable  or  shorter  that  the 
decoherence  time  of  the  medium.  The  new  numerical  ultrafast  pulse  propagation  solver  allows  us  to 
calculate  the  temporal-  and  spectral  shape  of  the  pulse  as  a  function  of  key  molecular  parameters,  such 
as  transition  frequency,  values  of  the  transition  dipole  moments  and  difference  of  the  permanent 
electric  dipole  moments  between  different  states,  as  well  as  the  rates  of  decoherence  and  energy 
relaxation  between  the  states.  Furthermore,  our  model  also  allows  taking  into  account  realistic 
resonance  line  shapes  of  both  one -photon  and  multi-photon  transitions.  We  demonstrate  the  utility  of 
our  numerical  algorithm  by  modeling  the  optical  power  limiting  performance  of  a  medium  comprising 
known  organic  chromophore  and  comparing  our  results  to  analytical  solutions  obtained  in  the  limit 
where  the  role  of  coherent  effects  is  minimal. 

Most  previous  works  addressing  the  modeling  of  pulse  propagation  in  multi-photon  optical  power 
limiting  medium  rely  on  the  circumstance  that  coherence  of  multi-photon  optical  transitions  is 
typically  rather  short-lived,  on  the  time  scale  of  100  fs  or  less.  If  the  duration  of  the  incident  pulses  is 
substantially  longer  than  the  decoherence  time,  then  it  is  possible  to  assume  that  the  2PA  process  takes 
place  practically  instantaneously.  This  so-called  “incoherent”  description  of  the  light-matter 
interaction  is  valid  if  the  following  two  key  conditions  are  satisfied:  (i)  the  rate  of  the  excitation  of  the 
2PA  transition  remains  relatively  low  compared  to  the  rates  of  decoherence  of  the  excited  states  and 
(ii)  the  bandwidth  of  the  pulse  is  narrow  compared  to  the  spectral  absorption  features  of  the 
chromophore,  including  one-photon  as  well  as  multi-photon  transitions.  In  this  case  the  light  pulse 
behaves  essentially  as  monochromatic  wave  and  the  excitation  processes  in  the  optical  limiting 
medium  may  be  adequately  described  in  terms  of  kinetic  equations  linking  the  population  in  different 
energy  eigenstates  with  the  rates  of  excitation  and  rates  of  energy  relaxation.  Due  to  relative  simplicity 
of  “incoherent”  optical  power  limiting,  this  approach  has  been  instrumental  in  analyzing  the 
performance  of  NLO  media  from  pico-  to  microsecond  time  scale,  including  those  which  take 
advantage  of  step-wise  absorption  from  excited  singlet-  and  triplet  states  [1-7]. 

In  the  project  described  here  we  are  addressing  what  we  call  “coherent”  regime  of  optical  limiting, 
where  the  assumption  of  instantaneous  transitions  may  no  longer  be  valid.  This  regime  is  relevant  in 
case  of  optical  limiting  media  comprising  organic  chromophores  will  distinct  electronic-vibrational 
spectral  features  and  if  the  incident  pulse  has  a  high  intensity,  Iin  >1-10  GW  cm'  and  the  duration  of 
the  pulse  is  100  fs  or  less.  In  this  case  the  ability  of  “incoherent”  models  to  predict  the  energy  of  the 
transmitted  pulse  is  rather  limited  because  kinetic  equations  do  not  adequately  describe  the  change  of 
the  temporal-  and  spectral  shape  of  the  pulse.  In  fact,  if  the  2PA  is  considered  instantaneous,  then  this 
is  equivalent  to  demanding  that  the  attenuation  of  the  intensity  at  any  time  is  strictly  proportional  to 
the  instantaneous  intensity  at  that  time.  This  means  that  the  phase  of  the  incident  pulse  has  no 
influence  on  the  optical  power  limiting,  which  is,  of  course,  because  the  kinetic  equations  account  only 
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for  the  intensity  and  not  the  phase.  In  order  to  take  adequately  into  account  coherent  effects  we  use  a 
more  accurate  photo-physical  model  that  is  based  on  solving  density  matrix  (DM)  equations  of  motion 
describing  how  the  incident  light  interacts  with  the  chromophores.  Even  though  the  solution  of  DM 
equations  is  computationally  much  more  demanding  than  solution  kinetic  equations,  this  approach 
allows  us  to  explicitly  accounts  for  the  temporal-  and  spectral  phase  of  the  pulse.  Previously,  Wang  et 
al.  [8]  have  solved  numerically  the  Maxwell-Bloch  DM  equations  using  finite-difference  time-domain 
(FDTD)  technique  to  study  the  propagation  of  femtosecond  laser  pulses  in  strong  two-photon 
absorption  organic  molecular  medium.  However,  most  likely  due  to  inherent  limitations  of  the  finite- 
difference  methods,  their  solutions  were  restricted  to  very  short  propagation  distances  less  than  10  qm 
and  extremely  short  pulses  containing  only  a  few  oscillation  cycles.  Our  calculation  uses  an  alternative 
approach  which  permits,  for  the  first  time,  the  study  of  the  spectral-,  temporal  and  phase  reshaping  for 
much  longer  pulses  containing  up  to  hundreds  of  optical  cycles  and  for  arbitrary  propagation  length. 

To  achieve  this  objective  we  have  combined  the  numerical  solution  of  DM  equations  with  a  modified 
step-by-step  description  of  the  propagation  in  the  forward  direction.  We  have  developed  a 
Mathematica- based  numerical  code  that  allows,  for  the  first  time,  accurate  quantitative  description  of 
the  temporal-  and  spectral  shape  of  the  pulse  as  it  propagates  in  the  2PA  medium  with  arbitrary  large 
optical  depth  and  arbitrary  long  propagation  distance.  Furthermore,  our  approach  allows  taking  into 
account  non-Forentzian  line  shape  of  the  resonance  transitions  (both  one-photon  as  well  as  multi¬ 
photon).  Because  the  Forentzian  line  shape  that  is  typically  used  in  modeling  2PA  has  far-reaching 
wings,  the  choice  of  actual  line  shape  may  strongly  influence  the  outcome  of  pulse  propagation 
especially  if  the  spectral  width  of  the  transitions  is  broad  on  the  order  of  1000  cm'1  or  more.  We 
address  this  issue  by  using  stochastic  density  matrix  (SDM)  model  that  has  been  recently  described  in 

[9] , 

The  remainder  of  this  report  is  organized  as  follows.  In  Section  3  we  present  theoretical  background, 
which  includes  a  simple  analytical  solution  for  plane  wave  optical  power  limiting  by  2PA,  description 
of  the  density  matrix  model  of  the  2PA  chromophores  and  description  of  the  the  step-by-step 
propagation  solution.  We  also  briefly  discuss  some  of  the  most  important  technical  aspects  of  the 
calculation  process.  In  Section  4  we  present  results  where  we  study  propagation  of  a  Gaussian-shaped 
incident  pulse  of  100-fs  duration  in  a  three-level  2PA  medium  modeled  after  a  well-known  2PA 
absorber.  We  carry  out  this  calculation  for  a  broad  range  of  input  peak  intensity  10-10  GW/cm  .  We 
use  the  calculation  data  to  evaluate  the  dependence  of  the  nonlinear  transmission  on  the  input  pulse 
energy  density  for  various  propagation  distances  in  the  medium.  Finally,  in  Section  5  we  present 
conclusions  and  discuss  outlook  for  future  applications. 

3.  Theoretical  background  and  implementation. 

3.1.  Simple  analytical  model  of  2PA-based  optical  power  limiting. 

Textbook  example  of  incoherent  2PA  consists  of  linearly  polarized  monochromatic  light  with  intensity 
that  is  constant  in  time.  The  change  of  the  intensity  with  the  propagation  distance,  dz,  may  be 
expressed  as: 

±Ll(z)  =  -\o2nl\z),  (1) 

dz  2 

where  02  is  the  molecular  2PA  cross  section  and  n  is  the  concentration  (number  of  2PA  chromophore 
molecules  per  unit  volume).  The  factor  Vi  is  introduced  as  a  matter  of  following  a  previous  convention 

[10] ,  The  differential  equation  leads  immediately  to  the  well-known  2PA  attenuation  law: 
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/(zWoO  +  ioiB/.z)-',  (2) 

where  Io  is  the  incident  intensity  (number  of  photons  per  second  per  cm  ).  As  an  illustration  of  equ. 
(2),  let  us  calculate  how  the  intensity  would  be  attenuated  if  we  assume  that  the  medium  comprises 
trans-4,4’-Bis(diphenylamino)stilbene  (BDPAS)  molecules  whose  experimentally  measured  2PA 
spectrum  [10]  is  shown  in  Figure  1. 


Laser  wavelength  ,  nm 

900  850  800  750  700  650  600 


Figure  1.  2PA  cross  section  spectrum  (dots)  of  trans-4,4’-Bis(diphenylamino)stilbene  (BDPAS)  dissolved  in  methylene 
chloride  [10].  The  spectral  shape  of  one-photon  absorption  transitions  plotted  vs  twice  the  normal  wavelength  is  shown  for 
comparison  (dashes  line).  Insert  shows  the  chemical  structure  of  BDPAS. 


Suppose  that  the  incident  wavelength  is  690  nm  corresponds  to  the  peak  2PA  cross  section  of  BDPAS, 
02  =  335  GM  (1  GM  =  IO'50  cm4)  and  that  the  chromophore  concentration  is,  n=0.5  1019  cm'3  (8.5 
mMol).  Left  panel  of  Figure  2  shows  the  relative  intensity  (normalized  to  the  incident  intensity)  as  a 
function  of  the  propagation  distance  z  (left  panel).  Right  panel  of  Figure  2  presents  the  same  data  but 
plotted  as  a  function  of  the  incident  intensity  for  different  z  values.  This  simple  model  predicts,  for 
eample,  that  for  the  maximum  propagation  distance  of  0.5  cm  the  output  intensity  will  be  clamped  at 
~70  GW  cm'2. 


Figure  2.  Attenuation  of  stationary  monochromatic  light  in  medium  comprising  BDPAS  molecules  at  concentration  n=0.5 
1019  cm"3.  Left  panel  -  intensity  transmission  as  a  function  of  the  propagation  distance  z  for  different  incident  intensity 
levels;  Right  panel  -  clamping  of  the  output  intensity  for  different  propagation  distances  0-0.5  cm. 
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3.2.  Density  matrix  equation  of  motion  for  three-level  system. 

In  our  approach  the  absorbing  medium  is  described  by  density  matrix,  p ,  and  the  electromagnetic 
field  is  given  by  real-valued  oscillating  electric  field  amplitude,  A(t,  z) .  The  change  of  the  state  of  the 
medium  with  time  is  obtained  by  solving  the  semi-classical  quantum-mechanical  DM  equation  of 
motion: 

d  A  1  /v 

—  pit)  =  —[H0  -  A(t, z)p, p(t)]  +  relaxation  terms,  (3) 

dt  i 

*  ^ 
where  H0  denotes  the  unperturbed  energy  of  the  molecule  and  p  is  the  electric  dipole  moment 

operator.  For  simplicity  we  are  assuming  that  the  optical  field  is  linearly  polarized  and  that  all  dipole 
moment  vectors  are  directed  in  the  same  direction  parallel  to  the  electric  field  vector.  At  this  stage  it  is 
sufficient  to  consider  that  the  chromophore  has  maximum  of  three  essential  energy  levels.  In  this  case 
the  three  operators  can  be  written  as  3x3  matrixes: 


(Poo 

Ai 

P02) 

'Poo 

Mi 

M2 N 

[Eo 

0 

0  "l 

A  o 

/V 

A  i 

Pn 

,  u  = 

Mo 

Mi 

M2 

and  H0  = 

0 

E\ 

0 

/V 

\P2  0 

Ai 

/V 

P22 ) 

Mi 

M2 ) 

1° 

0 

E>) 

of  which  first  two  are  complex  valued  and  last  is  real  valued.  Substituting  (4)  into  (3),  and  considering 
the  conservation  law,  pm(t)  +  pu(t)  +  p21(f)  =  f  gives  us  following  set  of  8  coupled  1st  order  ODE: 

A Rep0 ! (0  =  -}'() j Rep0 }  ( t )  +  [AU)Au() j  - Q V»'P() j (0  +  A(t)[Rep{ P>»P{)2 (0  +  ReP()2ll"Pl2(l>  ~  /m/'| 2ReP()2(l>  ~ Imp^Rep^U)  +  Imu() j (1  - 2 Pn<J)-p22 (/))] 
A Imp0 j (0  =  -y() !  Inip() j (t)  -[A(t)Ap0 }  {)\ReP{) ,  (0  +  AU)[-Imu}  pnp^O)  +  Imp^Imp^ 2 (!)  - Rep]  2Rep{)2 ( t )  +  Rep^Rep^ 2 (<)  -  Reu{) ^1-2 Pll(‘)~P22 (0)] 

A  Rep(p(t)  =  -y()2Rep(j2(t)  +  [A(t)Ai,()2 -0J2()]imp(j2U)  + A(t)[Rep]2Imp(jp)- Rep(jpnpn(t)  +  ImunRep()p)  - Imp()lRep]2U)  + - p^  pA-lp^U))] 
Almp()1(t)  =  -r02Imp02(t)-[A(t)Ap02  -  co^Rep^t )  +  AU)[Imu}  pmp{) p)-  } Impu(t)  -  Rep^Rep^p)  +  Reu^Rep^pt)  - Rep^Q.  - p}  p)  - 2 /=’22C0)] 

AReP]2(‘)  =  -YnRePu<‘)  +  M<')A/'|2  - co2l]Impu(t)  +  AU)[-Reu{pJmp{)p)  - Reu{)pnp{pU)  +  Imp(pRep()p)  +  Imp(pRep(pU)  +  Impn(pl  p)  - P22(t))] 
Almpl2(t)  =  -ruImpu(t)-[A(l)Apu  -  a>2l]Repu(t )  +  A(t)[-lnip(phnp()p)  +  lmp(plmp(p(t)  - Rep^Rep^p)  +  Rep^Rep^t)  -  Repn(px  p)  -  p22(t))] 
Aplp)  =  -Tl0plp)  +  r2lp22(t)-2A(t)[Rep0p)Imp0p)-Repu(t)Impu(t)-Imp0p)Rep0p)  +  Impu(t)Repu(t)] 

Ap22{t)  =  -(T2p  T2\PP22^  “  2A^Re^02^Imp02^  ~  Rep\ 2 ^ImP\ 2 (0  “ Im^2^Rep02^  +  Imp\2^Rep\2^ 

(5) 

Here  I\.  is  the  population  (energy)  relaxation  rate  and  y/;  (i  *  j)  is  the  rate  of  decoherence  between  the 

levels  i  and  j.  We  have  also  introduced  a  notation  for  the  transition  frequencies:  C0y  =  (E.  -Ej)/  h,  and 

for  the  dipole  moments:  A ptj  =  p2u  -  pu.  In  the  above  equations  the  numerical  value  of  the  Planck 

constant  is  absorbed  in  the  amplitude  A(t).  Typical  initial  conditions  for  the  DM  equations  is  when  at 
t=0,  i.e.  before  the  arrival  of  the  excitation  pulse  all  matrix  elements  other  than  the  ground  state 
population  are  equal  to  zero: 

^/>01(0)  =  ImpQl( 0)  =  Repu( 0)  =  Impl2( 0)  =  Rep^O)  =  Imp^O)  =  p{  ^0)  =  A»22 (0)  =  0; 

P00(0)  =  f 


(6) 
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Not  that  the  rotating  wave  approximation  (RWA)  commonly  used  in  solving  the  nonlinear-optical  DM 
equations  is  not  applicable  here.  This  is  because  accurate  description  of  2PA  and  other  essentially 
multi-photon  processes  depends  on  retaining  of  the  higher  harmonics  created  through  nonlinear 
interaction  of  the  oscillation  driving  field  with  the  molecular  system,  and  which  RWA  is  specifically 
designed  to  suppress  [11,  12].  For  a  closely  related  reason  we  cannot  use  the  common  slowly  varying 
envelope  approximation  (SWEA),  because  that  would  also  suppress  and  distort  2PA  in  the  coherent 
regime.  Nevertheless,  numerical  solution  of  (5)  (along  with  the  initial  conditions)  with  sufficiently 
high  accurate  is  quite  within  the  realm  of  capability  of  commercial  numerical  ODE  solvers,  especially 
if  the  applied  field  does  not  undergo  too  many  oscillation  periods.  As  will  be  discussed  in  more  detail 
below,  we  are  using  for  this  purpose  NDSolve  function  of  Mathematica  v.8  package. 


A  simple  example  of  the  excitation  pulse  is  that  with  carrier  frequency,  cop  =  2 jivp  =  2jic  /  Ap  and 
Gaussian  envelope  with  the  duration  xp  (intensity  FWHM): 


A(t)  =  Ape 


Tp  21n2 


cos  cop(t-tp), 


(7) 


where  Ap  is  the  pulse  peak  amplitude  and  tp  is  the  pulse  peak  delay.  Figure  3  shows  the  normalized 
amplitude  (left  panel)  and  the  power  spectrum  of  a  pulse  (right  panel)  with  the  parameter  values,  xp= 
100  fs  and  A p=  690  nm.  We  choose  these  parameters  because  for  low  peak  intensities  the  pulse 
spectrum  may  be  still  considered  as  nearly  monochromatic  (compared  to  the  width  of  the  2PA 
spectrum),  which  facilitated  comparison  between  numerical  solution  and  the  above  simple  analytical 
solution. 


Laser  wavdengfh  .  nm 

900  830  800  750  700  650  600 


400 


11000  12000  13000  14000  15000  16000  17000 

Freqierry  ,  cm-1 


Figure  3.  Normalized  time  domain  amplitude  (left  panel)  and  the  frequency  domain  power  spectrum  (right  panel)  of  a 
pulse  with  typical  the  parameter  values  used  in  the  calculation,  rp=  100  fs  and  Ap=  690  nm.  The  2PA  spectrum  of  BDPAS 
is  shown  for  comparison. 

3.3.  Simulation  of  realistic  non-Lorentzian  line  shape. 

The  phenomenological  relaxation  introduced  in  equ.  (5)  is  responsible  for  finite  spectral  line  widths  of 
the  transitions.  As  the  DM  is  written  now,  the  relaxation  is  due  to  the  constant  decoherence  rates,  Yu- 
Constant  relaxation  rate  that  does  not  depend  neither  on  time  or  frequency  correspond  to  Markovian 
process,  which  leads  directly  to  Lorentzian  line  shape.  Experiment  shows,  however,  that  actual  line 
shapes  lack  the  far-reaching  “wings”  that  are  characteristic  of  Lorentzian  function.  In  fact,  the  real  line 
shapes  such  as  those  shown  in  Figure  1  are  closer  to  Gaussian  (or  a  superposition  of  Gaussians).  This 
circumstance  is  particularly  critical  for  modeling  of  multi-photon  absorption  because  at  the  long 
wavelengths  the  Lorentzian  wing  of  the  1PA  resonance  may  still  be  strong  enough  to  interfere  with  the 
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2PA  transition,  which  leads  to  an  unrealistic  description  of  the  pulse  propagation  at  the  2PA 
resonance. 

We  have  previously  shown  that  realistic  non-Lorentzian  line  shapes  can  be  modeled  by  a  Monte  Carlo 
type  DM  calculation,  where  the  transition  frequency  contains  a  stochastic  component  with  certain 
amplitude  and  frequency  distribution  [9].  For  this  we  substitute  in  the  equ.  (5)  the  transition 
frequencies  between  energy  levels  i  and  j  with  the  following  expressions: 

a)ij^a).lj[\  +  doJij7-]..(t)l  (8) 

where  rjy(t)  is  the  stochastic  time-domain  “noise  function”  with  unit  maximum  amplitude  and  Scoy  is 
the  corresponding  amplitude.  The  power  spectrum  as  well  as  the  amplitude  of  the  stochastic  “noise 
function”  are  adjusted  to  achieve  best  correspondence  to  the  measured  1PA  and  2PA  spectra. 

Our  next  step  is  to  determine  molecular  parameters  for  the  three-level  model  that  best  describe  the 
experimental  properties  of  BDPAS.  The  corresponding  values  are  shown  in  Table  1.  Here  and  below 
we  assume  that  transition  dipole  moments  are  real  valued. 


Table  1. 
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24800  cm'1 

29050  cm'1 
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0.27 

Figure  4  shows  the  calculated  1PA  and  2PA  spectra  of  BDPAS  in  comparison  to  the  measured 
experimental  data.  Good  correspondence  is  achieved  on  both  linear  scale  (left  panel)  and  on 
logarithmic  scale  (right  panel).  The  quantitative  agreement  is  sufficiently  good  for  current 
demonstration  purposes  but  could  be  further  improved  by  fine  adjustment  of  the  parameters  of  the 
model. 


Figure  4.  Calculated  spectra  of  BDPAS  (orange)  compared  to  the  measured  experimental  data  (dark  blue);  Upper  panel  - 
1PA  spectrum;  Lower  panel  -  2PA  cross  section  spectrum.  The  parameter  values  used  are  shown  in  Table  1. 
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3.4.  Step-by-step  propagation  method. 

The  task  of  this  project  is  to  develop  numerical  solutions  that  describe  the  propagation  of  a  plane  wave 
front  linearly  polarized  ultrashort  light  pulse.  To  address  this  problem,  one  typically  would  needs  to 
solve  the  DM  equations  (5)  simultaneously  with  one-dimensional  Maxwell’s  equations: 

1  r)2  AjT  r)2 

W2E(t,z)--  —  E(t,z)  =  —  —  P(t,z).  (9) 

c  at  c  at 

Here  the  macroscopic  polarization  is  found  from  the  density  matrix: 

P(t,z)  =  n{Trjup(t,z)},  (10) 

where  the  averaging  <. .  .>  is  carried  out  by  integration  over  an  intrinsic  sub-wavelength  scale 
distribution  of  parameters.  Typical  distribution  function  g(co,Q)  describes  inhomogeneous  spectral 
shifts  of  transition  frequencies  and  different  orientations  of  the  molecule: 

(/(f,z»  -$g(co,Q)f(t,z)dcodQ .  (11) 

Unfortunately,  because  we  cannot  take  advantage  of  neither  RWA  nor  SWEA  approximations,  the 
exact  solution  constitutes  a  rather  formidable  computational  problem.  We  get  around  this  difficulty  by 
reducing  the  wave  equation  (9)  to  a  much  simpler  case,  where  reasonable  numerical  solution  is 
feasible.  We  note  that  in  complex  notation  the  amplitude  of  the  pulse  propagating  in  vacuum  in  the 
positive  z  direction  may  be  written  as  follows: 

&(t,z)  =  ^t-zlc)eico”(,-zlc\  (12) 

The  wavy  hat  (A=A%  in  Mac’s  rendering,  sorry)  over  the  amplitude  indicates  the  analytical  signal: 

00  00 

%%)  =  A(t)  +  illmj  ew!t  J  eioJt'  A(t')  dt'dco  ,  (13) 

0  -oo 

where  as  before  A(t)  stands  for  the  real  amplitude.  We  use  Fast  Fourier  Transform  (FFT)  and  inverse 
FFT  routines  in  the  Mathematica  package  to  implement  the  calculation  of  analytical  signals  equ  (13). 

Fet  the  pulse  be  incident  at  z=0  at  2PA  medium  consisting  of  chromophores  at  a  concentration  of  n 
molecules  per  unit  volume.  The  goal  is  to  find  the  amplitude,  |  M$,d)  ,  and  phase,  Arg[$t,d)\ ,  of  the 
pulse  at  the  propagation  distance  d. 

Fet  us  divide  the  medium  into  N  thin  slices  of  thickness  Az=d/N.  Each  slice  will  contain  only  a  small 
number  of  molecules  per  unit  area,  n  Az,  such  that  the  polarization  (scattering)  induced  in  each  slice 
will  have  only  a  small  effect  on  the  pulse  amplitude.  However,  with  the  propagation  through  many 
such  slices  the  effect  of  the  induced  polarization  accumulates,  which  in  turn  will  result  in  the 
attenuation  and  change  of  the  pulse  shape  (amplitude  and  phase).  This  scheme  is  illustrated  in  Figure 
5,  which  shows  radiation  pattern  of  dipoles  excited  in  the  nth  layer  by  the  incoming  plane  wave.  The 
dipole  emits  coherently  in  the  forward  direction.  Backward  emission  (and  backward  propagation)  such 
as  reflection  is  disregarded.  Only  coherent  emission  (scattering)  in  the  forward  direction  is  considered. 
The  amplitude  of  the  incident  wave  combines  coherently  with  the  (much  weaker)  forward  scattering 
amplitude.  The  resulting  coherent  superposition  of  the  two  amplitudes  creates  the  input  pulse 
amplitude,  which  serves  as  input  amplitude  for  the  (n+l)th  layer  etc. 
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plane  wave 


Figure  5.  Modified  split  step  method  for  the  pulse  propagation  in  2PA  medium.  Incident  plane  wave  travels  from  left  to 
right.  Induced  dipole  scatters  i.e.  emits  coherently  in  the  forward  direction.  The  scattering  consists  of  linear  component  and 
much  weaker  nonlinear  component,  where  only  the  last  is  related  to  2PA.  Combined  amplitude  of  the  incident  wave  and 
the  scattered  wave  constitute  the  input  wave  for  the  (n+l)th  layer. 

The  expression  for  the  (n+l)th  amplitude  is: 

(0  =  ^{t)  -ft  —  —  P\t) .  (14) 

N  dt 

The  real  polarization  is  calculated  as  follows: 

Pn(t)  =  n2[p^Rep^{t)  +  p^Rep^it)  +  ^Rep^t)],  (15) 

where  the  density  matrix  elements  are  found  by  solving  the  equ.  (5)  for  the  nth  layer.  The  constant  (5 
represents  other  factors,  such  as  index  of  refraction  of  the  medium,  orientation  of  the  dipole  moments 
etc.  which  are  not  taken  into  account  in  this  model.  The  value  of  /3  is  determined  in  an  empirical 
manner  e.g.  by  taking  the  pulse  carrier  frequency  at  the  peak  of  the  1PA  absorption  and  then  matching 
the  attenuation  at  low  incident  pulse  intensity  with  the  known  linear  absorption  coefficient  and  at 
given  chromophore  concentration.  In  addition,  in  our  model  we  apply  spectral  filtering  of  the 
polarization  amplitude  and  temporal  filtering  of  the  pulse  amplitude.  The  spectral  filtering  is  used  to 
eliminate  non-physical  effects  such  as  second  harmonic  generation.  The  temporal  filtering  is  used  to 
cuts  off  computational  artifacts  at  the  start  edge  and  the  end  edge  of  the  calculated  time  interval,  and 
which  occur  because  the  length  of  the  time  interval  is  finite. 

The  accuracy  of  the  step-by-step  propagation  method  can  be  increased  by  increasing  N.  However,  the 
nonlinear  polarization  that  is  responsible  for  the  2PA  effect  constitutes  only  a  small  fraction  of  the 
overall  induced  oscillating  dipole.  Consequently,  even  at  high  intensity,  most  part  of  the  polarization  is 
still  due  to  the  one  photon  resonance,  even  though  the  pulse  frequency  may  be  detuned  far  from  the 
1PA  peaks.  This  means  that  in  reality  to  achieve  reasonable  accurately  for  a  practical  level  of  2PA 
optical  limiting  (e.g.  50%  or  less  transmission  decrease)  the  number  of  the  slices  N  should  be  quite 
large,  which  in  turn  increases  proportionally  the  calculation  time  required.  We  have  addressed  this 
issue  by  solving  the  DM  equation  twice  for  each  layer:  one  time  with  normal  parameters  and  the 
second  time  with  the  2PA  switched  off  such  that  only  the  1PA  resonances  are  present.  The  last  is 
achieved  by  setting  to  zero  the  molecular  parameters  which  are  specifically  responsible  for  the  2PA. 
For  example,  if  we  set  p2i=0  and  Auio=0,  then  2PA  vanishes  while  the  “linear”  polarization  due  to  the 
0  — ►  1  resonance  transition  is  still  present.  We  subs  tract  the  “linear”  polarization  from  the  total  (linear 
+  nonlinear)  polarization,  which  gives  us  the  desired  2PA  effect  but  without  the  large  one-photon 
background.  The  described  method  allowed  us  to  substantially  increase  the  maximum  absorption  per 
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one  layer  while  maintaining  sufficient  accuracy  for  evaluation  of  optical  limiting.  This  in  turn  allowed 
us  to  decrease  the  overall  number  of  layers  needed  N  at  least  by  a  factor  of  10  and  accordingly 
decrease  the  calculation  time. 

Another  important  parameter  that  factors  into  the  overall  calculation  time  is  the  number  of  averages 
for  stochastic  DM  equation.  The  simulated  spectra  (without  propagation)  presented  in  Fig.  4  consists 
of  120  spectral  points,  where  each  represent  an  average  over  280  calculations.  The  total  time  for  this 
calculation  was  about  270  minutes  (~4.5  hours)  on  a  14  core  PC,  which  means  that  the  time  needed  to 
solve  the  stochastic  DM  for  one  pulse  carrier  frequency  is  about  2  min  (up  to  10  min  if  high  incident 
intensity.  If  we  consider  that  typical  number  of  layers  is  N=  1000  -  2000  and  that  solving  DM  for  each 
layer  requires  at  least  2  minutes,  then  the  overall  duration  of  the  calculation  would  be  still 
prohibitively  lengthy.  Fortunately,  we  found  a  way  how  to  shorten  the  calculation  further.  To  our 
pleasant  surprise  we  found  that  the  stochastic  averaging  works  equally  well  not  only  if  performed 
separately  for  each  layer,  but  also  if  it  is  performed  in  the  “time  forward”  direction.  From  the  point  of 
view  of  smoothing  out  the  stochastic  fluctuations  in  the  pulse  shape,  propagating  the  pulse  through  10 
layers  has  a  similar  effect  as  performing  the  averaging  over  10  samples  in  one  layer.  The  critical 
condition  is  of  course  that  the  effect  of  each  layer  remains  relatively  small.  With  this  in  mind,  we 
found  that  optimum  results  in  terms  of  reducing  the  calculation  time  versus  obtaining  a  sufficiently 
smooth  pulse  shape  are  achieved  if  the  averaging  in  each  layer  is  performed  about  10-20  time  (as 
opposed  to  280  times  in  Fig.  4).  This  gave  us  another  order  of  magnitude  reduction  in  the  speed  of  the 
calculation. 

4.  Results. 

4.1.  Verification  of  the  numerical  solver  by  comparison  to  known  analytical  solution. 

The  first  objective  is  to  show  that  at  relatively  low  input  intensity  the  nonlinear  transmission 
corresponds  to  the  simple  analytical  formula  equ.(2).  This  comparison  will  allow  us  to  verify  if  our 
calculation  is  indeed  delivering  qualitatively  and  quantitatively  reliable  results.  Figure  6  shows  the 
normalized  intensity  as  a  function  of  the  propagation  distance  for  the  incident  intensity  values  47  GW 
cm'  and  94  GW  cm'  .  The  concentration  of  BDPAS  molecules  is  n=0.5  10  cm'  in  all  cases.  At  these 
moderate  intensities  the  contribution  of  saturation  should  be  minimal.  Because  the  spectral  width  of 
the  100-fs  pulse  is  several  times  less  than  the  line  width  of  the  2PA  transition,  the  contribution  of 
coherent  effects  should  be  also  small.  Solid  line  is  the  numerical  result  and  the  dashed  line  is  the 
analytical  solution.  The  maximum  propagation  length  is  1  cm,  which  is  divided  into  N=2000  layers. 


z,cm 


Figure  6.  Normalized  transmitted  intensity  calculated  numerically  (solid  line)  and  analytically  (dashed  line)  plotted  as 
function  of  the  propagation  distance  for  two  incident  peak  intensities  47  GW  cm"2  and  94  GW  cm'2.  The  medium  consists 
of  BDPAS  molecules  at  concentration  n=0.5  1019  cm"3.  The  number  of  layers  N=2000  corresponds  to  the  maximum  path 
length  1  cm.  Right  panel  zooms  in  on  the  short  distance  section  0<z<0.05cm. 
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The  fact  that  our  numeric  result  corresponds  well  to  the  analytical  result  confirms  that  at  least  in  the 
low  to  moderate  intensity  regime  our  solver  is  reliable.  The  slightly  higher  transmission  for  the 
numerical  value  is  most  likely  because  analytical  formula  assumes  constant  intensity.  Note  that  the 
numerical  curves  show  slight  random  oscillations  or  “noise”  which  is  due  to  the  stochastic  fluctuations 
built  into  the  DM  equations.  It  is  an  interesting  observation  that  the  “noise”  does  not  increase  with 
larger  z  but  rather  tends  to  average  out.  This  property  is  in  fact  very  useful  because,  as  discussed  in  the 
previous  section,  this  allows  us  to  reduce  the  number  of  averaging  steps  per  each  layer.  This  in  turn 
provides  substantial  economy  in  terms  of  overall  computation  time. 

4.2.  Optical  clamping  and  saturation. 

As  the  next  step  we  increase  the  intensity  of  the  incident  pulse,  while  keeping  all  other  parameters  the 
same.  In  our  case  the  optical  power  limiting  is  accompanied  by  the  change  of  the  pulse  shape  (see 
below).  However,  to  facilitate  direct  comparison  to  the  analytical  model  shown  in  Figure  2,  we 
assume  here  that  pulse  shape  does  not  change,  and  scale  our  data  considering  Gaussian  shape  (for  100- 
fs  duration  Gaussian  pulse  1  GW/cm  peak  intensity  corresponds  to  energy  density  ~  106.4  uJ  cm' ). 
Figure  7  shows  the  effective  pulse  intensity  as  a  function  of  the  incident  peak  intensity  for  different 
propagation  propagation  distances,  z=  0  -  0.5  cm.  The  data  shows  that  at  given  concentration  n=0.5 
109  cm'3  the  minimum  propagation  distance  is,  z>  0.4  -0.5  cm.  The  effective  clamping  intensity  is 
about  100  GW  cm'  ,  which  compares  well  to  the  value  70  GW  cm'  obtained  above  for  the  cw  case 
(dashed  curve  in  Fig.  7).  On  the  other  hand,  if  the  incident  intensity  is  higher  than  >10  GW  cm'  ,  then 
the  transmission  increases  again  due  to  the  saturation  of  the  absorption. 


Figure  7.  Intensity  of  100-fs  Gaussian  pulse  propagating  in  the  2PA  medium  (BDPAS  at  concentration  n=0.5  1019  cm"3). 
Left  panel  -  effective  pulse  peak  intensity  plotted  as  the  function  of  the  incident  intensity  for  different  propagation 
distances,  z  =  0  -  0.5  cm.  Right  panel  zooms  in  on  the  larger  propagation  distances  showing  optical  clamping.  Dashed 
curve  -  data  according  to  the  analytical  solution  for  cw  light  and  z=0.5  cm  (see  Figure  2). 


We  should  mention  that  at  parameter  values  used  here  the  saturation  effects  may  be  also  quantitatively 
described  by  incoherent  rate  equations  [3,6,7].  We  are  in  the  process  of  performing  such  comparison 
calculations  to  further  verify  the  validity  of  our  new  solver.  In  conclusion  of  the  above  section  we  may 
say  that  we  have  verified  that  our  numerical  method  indeed  works  and  is  able  to  quantitatively  model 
the  performance  of  realistic  optical  power  limiting  media. 

4.3.  Coherent  propagation  effects.  Time-  and  spectral  domain  pulse  shaping. 

Now  we  can  finally  come  to  the  results  that  show  how  the  propagation  changes  the  temporal  shape,  the 
power  spectrum  of  the  pulse.  We  underline  that  the  results  presented  below  are  new  and  can  be 
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obtained  only  with  full  DM  solution,  and  cannot  be  obtained  neither  with  analytical  nor  with  less 
sophisticated  numerical  models.  Figure  8  shows  the  change  of  the  shape  of  the  pulse  as  a  function  of 
propagation  distance  for  two  different  input  peak  intensity  values.  In  the  case  of  moderate  intensity  (47 
GW  cm  )  the  time  domain  peak  diminishes,  while  the  overall  shape  remains  symmetrical.  This 
corresponds  to  lengthening  of  the  pulse  (sometimes  called  flat-topping).  In  the  spectral  domain  the 
shape  remains  more  or  less  constant;  No  new  frequency  components  are  generated.  At  higher  intensity 
(250  GW  cm  )  substantial  change  of  the  temporal  and  spectral  shape  occurs.  In  the  time  domain,  the 
peak  shifts  to  later  times.  Most  significantly,  the  peak  of  the  power  spectrum  shifts  to  higher 
frequency.  At  even  higher  intensity  the  spectrum  develops  oscillation  side  components. 


-200.  -100.  0.  100.  200. 


t,  fs 


-200.  -100.  0.  100.  200. 


t,fs 


Figure  8.  Change  of  pulse  shape  as  function  of  propagation  distance  0  -0.25  cm  for  input  peak  intensity  47  GW  cm2  (left 
panel)  and  250  GW  cm2  (right  panel).  BDPAS  concentration  n  =  0.5  1019  cm"3.  Maximum  propagation  distance  z=0.25  cm. 

Such  characteristic  drastic  changes  of  the  pulse  spectrum  propagating  in  the  2PA  medium  may  have  in 
fact  been  observed  experimentally.  Figure  9  compares  our  calculated  data  in  BDPAS  at  concentration 
n  =  0.5  10  cm'  with  the  measurements  performed  with  60-fs  pulses  at  675  nm  (estimated  peak 
intensity  4.2  104  GW  cm'3)  below-band-gap  propagating  through  CdS  quantum-dot-doped  waveguides 

[13]. 
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Figure  9.  Comparison  between  calculation  and  earlier  experiments  by  Guerreiro  et  al.  [13],  The  pulse  shapes  before 
(dashed)  and  after  (solid)  propagation  through  2PA  medium.  Left  panel  -  Time  domain  pulse  cross  correlation  (upper 
panel)  and  power  spectrum  (lower  panel).  The  experiments  were  performed  in  CdS  quantum-dot-doped  waveguides  with 
60-fs  pulses  at  estimated  peak  intensity  4.2  104  GW  cm'3. 

4.4.  Nonlinear  phase  change. 

Our  calculation  allows  determining  how  the  phase  of  the  pulse  changes  with  the  intensity  and  the 
propagation  distance.  This  information  can  be  used  to  study  the  effective  nonlinear  index  of  refraction 
caused  by  the  2PA  resonance  [14,  15].  Here  we  mention  only  one  illustration  of  this  possibility  leaving 
more  detailed  study  for  the  future.  Figure  10  shows  the  frequency  domain  phase  of  10-fs  pulse  after 
propagating  through  0.5  cm  of  BDPAS  medium  at  chromophore  concentration  n=0.5  10  cm'  for 
four  different  input  peak  intensity  values. 


Frequency,  cm  1 


Frequency,  cm 


Frequency,  cm 


Figure  10.  Frequency-domain  phase  of  10-fs  pulse  after  propagating  through  0.5  cm  of  BDPAS  medium  at  chromophore 
concentration  n=0.5  1019  cm"3  for  input  peak  intensity  values,  47,  188,  752  and  3008  GW  cm'2.  Orange  dots  -  phase  change. 
Dashed  line  -  spectral  amplitude  of  the  incident  pulse.  Solid  line  -  the  spectral  amplitude  after  propagation. 
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At  moderate  intensity  (upper  left  panel)  the  phase  change  is  relatively  small  and  close  to  constant  for 
all  frequency  components  of  the  pulse.  When  the  intensity  increases,  so  does  the  phase  change.  At  the 
higher  intensity  (two  lower  panels)  the  spectral  phase  develops  a  discontinuity  (jt-jump).  The 
frequency  of  the  phase  jump  coincides  with  the  frequency  where  the  spectral  amplitude  changes  sign 
and  intensity  goes  to  zero.  This  behavior  is  reminiscent  of  the  phase  change  of  soliton  amplitude  that 
develop  when  a  coherent  pulse  propagates  in  one  photon  resonance  medium.  Possibility  of  formation 
of  solitons  in  the  2PA  resonance  medium  has  been  considered  previously  in  [16-18]. 

5.  Conclusions  and  outlook. 

We  have  developed  a  versatile  numerical  technique  that  allows  us  to  calculate  the  energy,  the  temporal 
shape,  the  spectrum  and  the  phase  of  an  ultrashort  pulse  propagating  in  multi-photon  absorbing 
medium  with  fully  realistic  spectral  shapes  and  absorption  profiles.  Our  approach  is  based  on  solving 
straightforward  ODE  instead  of  Maxwell-Bloch  type  equations,  which  usually  require  FTDT 
calculations.  Furthermore,  our  solver  is  implemented  using  100%  commercial  software  ( Mathematica 
ver.  8)  running  on  a  conventional  multi-core  PC  platform,  which  avoids  complications  and  increased 
costs  associated  with  need  for  custom-produced  code  and/or  large-scale  computing  facilities.  This  fact 
in  itself  is  quite  remarkable,  especially  because  we  are  solving  a  set  of  stochastic  equations  that  require 
extensive  averaging.  First  we  verify  the  accuracy  of  our  calculation  by  comparing  to  some  simple 
analytical  formulas  and  then  we  demonstrate  its  utility  by  analyzing  the  optical  power  limiting 
performance  of  a  known  two-photon  organic  chromophore.  We  also  show  that  the  calculation  yields 
full  information  about  the  amplitude  and  phase  of  the  pulse. 

booking  forward,  we  envisage  that  our  new  solver  can  be  applied  for  the  simulation  and  analysis  of  a 
broad  range  of  physical  problems  including  coherent-  and  incoherent  regimes  of  optical  power 
limiting,  saturation,  CEP  effects,  soliton  formation  etc.  It  can  be  also  used  for  optimization  of 
materials  that  are  used  in  multi-photon  microscopy  and  imaging. 
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7.  List  of  Symbols,  Abbreviations  and  Acronyms 

02  two  photon  absorption  cross  section 

Ujj  electric  dipole  moment  of  transition  from  state  i  to  state  j 

Ajiij  permanet  electric  dipole  moment  difference  between  state  i  and  state  j 

2PA  two  photon  absorption 

BDPAS  trans-4, 4  ’-Bis(diphenylamino)stilbene 

SWEA  slowly  varying  envelope  approximation 

RWA  rotating  wave  approximation 

ODE  ordinary  differential  equations 

FDTD  finite  difference  time  domain 

DM  density  matrix 

SDM  stochastic  density  matrix 

CEP  carrier-envelope  phase 
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8.  Appendix 

Mathematica  program  listing. 


(★  Calculate  absorption  and  shape  of  a  pulse  propagating  in  three- level  2 PA  medium  *) 

(*  Aleks  Rebane  April  2011.  All  rights  reserved  *) 

Off [General: :"spelll"] 

(*  Set  directory  to  write  data  files  *) 

ResetDirectory [] ; 

SetDirectory [uRebane/Math/Y2011/Y20110614" ] ; 

Directory [] ; 

TableForm [FileNames [ ] ] ; 

(*  Check  for  available  kernels  *) 

kernelsOO  =  ParallelEvaluate [$KernelID] 

kernelsOl  =  {1,  2,  3,  4,  5,  6,  7 ,  8,  9,  10,  11,  12,  13,  14,  15,  16,  17,  18,  19,  20,  21,  22,  23,  24}; 
kernels02  =  Dimensions [kernelsOl]  [ [1] ] ; 

(*  01.  Pulse  peak  amplitude  20000  esu;  Write  data  to  "data_01x.txt"  *) 

(*  Master  kernel  section  *) 

(*  Define  initial  distribured  functions  and  parameters  *) 

noisefilterOO  =  Table  [Exp  [- 0 . 0032  *  (k  -  1)  A2]  +  Exp[-0.0032*  (k  -  1202)  A2]  ,  {k,  -1,  1202,  1}]  ; 
equOO  :=  { 

repOl'  [t]  ==  -yOl  *  repOl  [t]  +  (reAO  ApOl  -  o)10)  impOl  [t]  + 

reAO  (repl2  imp02  [t]  +  rep02  impl2  [t]  -  i_mpl2  rep02  [t]  -  imp02  repl2  [t]  -  impOl  pll  [t]  + 
impOl  (1  -  pll  [t]  -  p22  [t] )  )  , 

impOl' [t]  ==  -yOl  *  impOl  [t]  +  (- reAO  ApOl  +  o)10)  repOl  [t]  + 

reAO  (- impl2  imp02  [t]  +  imp02  i_mpl2  [ t]  -  repl2  rep02  [ t]  +  rep02  repl2  [ t]  +  repOl  pll  [t]  - 
repOl  (1  -  pll  [t]  -  p22  [t] ) )  , 
rep02'[t]  ==  -y02  *  rep02  [t]  +  (reAO  Ap02  -  <i)20)  imp02[t]  + 

reAO  (repl2  impOl  [t]  -  repOl  impl2  [t]  +  impl2  repOl  [t]  -  impOl  repl2  [t]  +  im^02  (1  -  pll  [t]  -  p22  [t] )  - 
imp02p22[t] )  , 

imp02'[t]  ==  -y02  *  imp02  [t]  +  (- reAO  A/i02  +  w20)  rep02  [t]  + 

reAO  (impl2  impOl  [t]  -  impOl  impl2  [t]  -  repl2  repOl  [t]  +  repOl  repl2  [t]  -  rep02  (1  -  pll  [t]  -  p22  [t] )  + 
rep02p22[t] )  , 

repl2'[t]  == -yl2  ★  repl2  [t]  +  (reAO  Apl2  -  w21)  impl2[t]  + 

reAO  (-rep02  impOl  [t]  -  repOl  imp02  [t]  +  imp02  repOl  [t]  +  impOl  rep02  [t]  +  impl2  pll  [t]  -  impl2  p22  [t] )  , 
impl2'[t]  ==  -yl2  ★  impl2  [t]  +  (- reAO  Apl2  +  w21)  repl2[t]  + 

reAO  (-imp02  impOl  [t]  +  impOl  imp02  [t]  -  rep02  repOl  [t]  +  repOl  rep02  [t]  -  repl2  pll  [t]  +  repl2  p22  [t] )  , 
pll'  [t]  =  -yll  *  pll  [t]  +  y22  *  p22  [t]  +  2  reAO  (-repOl  impOl  [t]  +  repl2  impl2  [t]  +  impOl  repOl  [t]  -  impl2  repl2  [t] )  , 
p22'  [t]  ==  -y22  *  p22  [t]  +  2  reAO  (-rep02  imp02  [t]  -  repl2  impl2  [t]  +  imp02  rep02  [t]  +  impl2  repl2  [ t ]  )  , 
repOl  [0]  =  0,  impOl  [0]  =  0,  rep02  [0]  =  0,  imp02  [0]  =  0,  repl2  [0]  ==  0,  impl2  [0]  =  0,  pll  [0]  ==  0,  p22  [0]  ==0}  /  . 
o)21  -►  (o)20  -  ojIO)  ; 

AO  =  ampOl  *  Exp  [-  ( (t  -  tl)  *  Sqrt  [2  .  *  Log  [2]  ]  /  dt)  A  2]  ★  Cos  [w  *  (t  -  tl)  ]  ; 


paraOO  =  | 

ampOl ->  20  000  .  ,  (*  pulse  peak  amplitude  ESU  *) 
tl ->  300 .  ,  (*  pulse  peak  time  fs  *) 

dt->100.,  (*  pulse  duration  fs  *) 

tranfreqlO ->  24  800  .  x  2  7r  *  3  .  x  10"5  ,  (*  0->l  frequency  cm-1 

tranfreq20 ->  29  050  .  x  2  7r  *  3  .  x  10-5  ,  (*  0->2  frequency  cm-1 

w ->  14  525  x  2  7r  *  3  .  x  10-5  ,  (*  pulse  carrier  frequency  cm-1 

<5o)10-»0.2,  (*  noise  amplitude  0-*l  transition  ★) 
<5o)20-»0.27  (*  noise  amplitude  0->2  transition  *) 


*) 

*) 

*) 
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paraOl =  { 

(*  Transition  dipole  moments,  Debye  *) 

re/vOl  -4  (6.0)  *  2  tt  *  1 . 509  x  10  A  -  7  , 
im^  0 1  -4  0  .  , 

re*/12-*  (9.0)  *  2  tt  *  1 . 509  x  10  A  -  7  , 
im/il2  -4  0  .  , 

reA/02  -4  0 , 

im£/02  -4  0  .  , 

(*  Permanent  dipole  moments,  Debye  *) 

A/./ 01  -4  (3.4)  *  2  7r  *  1 .509  x  10A -7, 

A/il2  -4  0.0, 

Aa/02  -4  0  .  , 

(*  Relaxation  rates,  fs-1  *) 

yOl  -4  2  tt*  0.00001, 
yl2->  2  7r  *  0 . 0001 , 
y02-*  2  tt*  0.0001, 
yll  -4  2  7r  ★  0 . 000001 , 
y22  ->  2tt  *  0.001 
} ; 

Al  =  A0  /  .  paraOO; 
equOl  =  equOO  /.  paraOl; 

DistributeDefinitions [noisefilterOO ,  equOO,  Al,  paraOO,  paraOl,  para03] ; 

(*  Parallel  evaluation  starts  using  "kernelsOl"  *) 
timeOO  =  ParallelEvaluate [TimeUsed [ ] ,  kernelsOl] ; 

(*  Create  blank  data  files  *) 
energydataO 1 = Table [ 0 ,  {k,  1,  1000}]; 
pulsedataOl  =  Table [0,  {k,  1,  10}] ; 
pulsedata02  =  Table [0,  (k,  1,  10}]; 
pulsedata03  =  Table [0,  {k,  1,  10}]; 
pulsedata04  =  Table [0,  (k,  1,  10}]  ; 
pulsedata05  =  Table [0,  (k,  1,  10}]  ; 

(*  Start  k  loop;  Propagate  k  layers  *) 

For[k  =  0,  k  <  1000,  k+4; 

DistributeDefinitions [Al] ; 

ParallelEvaluate [ 

(*  Create  dummy  data  files  *) 
da taO Oa = Table [0,  {t,  0,  1200,  1}]; 
dataOOb  =  Table [0,  {t,  0,  1200,  1}]; 
da taO Oc = Table [0,  {t,  0,  1200,  0.1}]; 

(*  Start  j  loop;  Solve  ODE  j  times  *) 

For [j  =  l,  j<2,  j+4; 

noiseO  la  =  Table  [RandomReal  [NormalDistribution  [0 ,  1.0]],  {k,  -1,  1202,  1}]; 
noise02a  =  Re [InverseFourier [noisefilterOO  ★  Fourier [noiseOla] ] ] ; 
noisefunca =  Interpolation [Table [ {k  -  1  ,  noise02a [ [k] ] } ,  [k,  1,  1204,  1 }  ]  ]  ; 
para02a  =  {&>10  tranfreqlO  *  (  1  +<5wl0  ★  noisefimca  [ t]  )  }  ; 

noiseOlb  =  Table  [RandomReal  [NormalDistribution  [0 ,  1.0]],  {k,  -1,  1202,  1}]; 
noise02b  =  Re [InverseFourier [noisefilterOO  *  Fourier [noiseOlb] ]]; 
noisefuncb =  Interpolation [Table [ {k  -  1  ,  noise02b [ [k] ] } ,  {k,  1,  1204,  1 }  ]  ]  ; 
para02b  =  {d)20  ->  tranfreq20  *  (  1  +6u)20  ★  noisefuncb  [t]  )  }  ; 
equOl  =  equOO  /  .  para02a  /  .  para02b  /  .  paraOO; 

(*  Solve  ODE  with  2PA  *) 

equ02  =  equOl  /  .  paraOO  /  .  paraOl  /  .  reAO  -►  Al  ; 

solOla =  NDSolve [equ02 ,  {rep01[t],  imp01[t],  rep02[t],  imp02[t],  repl2[t],  impl2 [ t] ,  pll [t] ,  p22 [t] } , 
{t,  0,  1200}  ,  MaxSteps  ->  500  000,  AccuracyGoal  ->  10,  PrecisionGoal  -4  10]  [  [1]  ]  ; 

(*  Repeat  the  solution  of  the  same  ODE  with  2PA  set  to  zero  ★) 

equOl  =  equOO  /.  para02a  /.  para02b; 

equ02  =  equOl  /  .  paraOO  /  .  para03  /  .  reAO  -4  Al  ; 

solOlb =  NDSolve [equ02,  {rep01[t],  imp01[t],  rep02[t],  imp02[t],  repl2[t],  impl2 [t] ,  pll [ t] ,  p22 [t] } , 
{t,  0,  1200}  ,  MaxSteps  -4  500  000,  AccuracyGoal  -►  10,  PrecisionGoal  -»10][[1]]; 

(*  Add  up  solutions  from  each  run  j  *) 

(*  Population  in  the  excited  states  1  and  2  *) 

datapopll  =  (pll[t]  /.  solOla)  -  (pll[t]  /.  solOlb) ; 
datapop22  =  (p22[t]  /.  solOla)  -  (p22[t]  /.  solOlb); 
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(*  Total  polarization  *) 
datapolOl = 

(2  (im*/01  *impd[t]  +  im^02  *imp02[t]  +  im^l2  *impl2[t]  +  re^Ol  * repOl [t]  +  re^02  *rep02[t]  + 
re*/12  *repl2[t])  /.  paraOl  /.  solOla)  - 

(2  (imA/01  *imp01[t]  + im^02  *imp02[t]  +  im^l2  *impl2[t]  +  re^Ol  *  repd  [t]  +  re^02  *rep02[t]  + 
re/il2  *repl2[t])  /.  para03  / .  solOlb); 
dataOla  = Table [da tapopll ,  {t,  0,  1200,  1}  ]  ; 
dataO lb  =  Table [datapop22,  {t,  0,  1200,  1}]; 
dataOlc = Table [datapolOl,  {t,  0,  1200,  0.1}]; 

(*  Add  results  of  j  runs  *) 
dataO  0 a  =  dataO  0 a  +  da  tad  a; 

(*  dataO Ob = dataO Ob + dataO lb;  *) 
dataOOc  =  dataOOc  +  dataOlc 
], 

(*  End  j  loop  *) 

kernelsOl  ]  ; 

(*  End  parallel  ODE  *) 

timed  =  ParallelEvaluate [TimeUsed [ ] ,  kernelsOl] ; 

Print [k] ; 

(*  Collect  data  and  average  *) 

outputOla  =  Total [ParallelEvaluate [dataO Oa,  kernelsOl] ]  /  4  .  ; 
outputdb  =  Total [ParallelEvaluate [dataO Ob,  kernelsOl] ]  /  4  .  ; 
outputOlc  =  Total [ParallelEvaluate [dataOOc,  kernelsOl] ]  /  4  .  ; 

(*  Evaluate  Re  polarization  *) 
repol  =  outputOlc; 

(*  Remove  polarization  spectral  components  outside  of  the  pulse  frequency  spectrum  *) 

complexpolspec  =  Inver seFourier [repol] ; 

Do [complexpolspec [ [i] ]  =0,  {i,  600,  12  001}]; 

Do [complexpolspec [ [i] ]  =0,  {i,  1,  300}]; 

(*  Cut  off  time-domain  artifacts  of  the  complex  polarization  before  and  after  the  pulse  *) 
list003  =  Fourier [complexpolspec] ; 

Do [list003 [ [i] ]  =0,  {i,  11500,  12001}]; 

Do [list003 [ [i] ]  =0,  {i,  1,  500}]; 

func003  =  Table [{0.1  (i  -  1) ,  Re [list003 [ [i] ] ] } ,  {i,  1,  12  001}] ; 
func004  =  Interpolation [func003]  [t] ; 
func005  =  D[func004,  t]  ; 

complexpolOl = Table [func005,  {t,  0,  1200,  0.1}]; 
complexpolspecOl = InverseFourier [complexpolOl] ; 

(*  Calculate  new  pulse  amplitude  *) 
oldpulsed  =  Table f Al .  It.  0.  1200,  0.1}  1  ; 


oldpulsespecOl  =  InverseFourier  [oldpulsed]  ; 

newpulsespec02  =  2  *  oldpulsespecOl  -  1 .  x 10 A  9  * complexpolspecOl; 
newpulsespec03  =  newpulsespec02; 

Do [newpulsespec03 [ [i] ]  =0,  {i,  600,  12001}]; 

Do [newpulsespec03 [ [i] ]  =0,  {i,  1,  300}]; 
newpulse03  =  Fourier [newpulsespec03] ; 

Do [newpulse03 [ [i] ]  =0,  {i,  1,  500}]; 

Do [newpulse03 [ [ i ]  ]  =0,  {i,  11500,  12  001}]; 

newpulse04  =  Table [{0.1  (i  -  1) ,  Re [newpulse03 [ [i] ]  ]  }  ,  {i,  1,  12  001}] ; 

Al  =  Interpolation [newpulse04]  [t] ; 

(*  Save  pulse  envelope  data  *) 

(*  Print [k] ;  *) 

energydataOl [ [k] ]  =  Total [Abs [newpulse03] A  2]  10  A - 11 ; 

If  [ 

Element [ (k  +  99)  / 100,  Integers] , 

1  =  (k  +  99)  /  100; 

pulsedataOl [ [1] ]  =  newpulse03; 

pulsedata02 [ [1] ]  =  newpulsespec03 [ [Range [300,  600] ] ] ; 
pulsedata03[ [1] ]  =  outputOla; 
pulsedata04  [  [1]  ]  =  outputdb; 
pulsedata05 [ [1] ]  =  complexpolOl 
,  0] 

(*  Print [Total [Abs [newpulse03] A2] ] ;  *) 

]  ;  (*  End  k  loop  *) 

time02  =  ParallelEvaluate [TimeUsed [ ] ,  kernelsOl] ; 

energydataOl  >>  dataO  la .  txt;  (*  Sum  of  energy  of  the  new  pulse  *) 
pulsedataOl  >>  da tadb .  txt;  (*  complex  time  amplitude  of  the  new  pulse  *) 
pulsedata02  >>  dataOlc . txt;  (*  complex  spectral  amplitude  of  the  new  pulse  *) 
pulsedata03  >>  datadd.  txt ;  (*  dynamics  of  pll  population  in  time  *) 
pulsedata04  >>  datade .  txt ;  (*  dynamics  of  p22  population  in  time  *) 

pulsedata05  >>  data01f.txt;  (*  Time  derivative  of  spectrally-  and  time  filtered  real  polarization  *) 

TableForm  [  { timeOl  -  timeOO,  time02  -  timeOl}  /  60] 


